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We present the results of linear analysis and two-dimensional local magneto- 
hydrodynamic (MHD) simulations of the Parker instability, including the effects 



of cosmic rays (CRs), in magnetized gas disks (galactic disks). As an unper- 
turbed state for both the linear analysis and the MHD simulations, we adopted 
an equilibrium model of a magnetized two-temperature layered disk with constant 
gravitational acceleration parallel to the normal of the disk. The disk comprises a 
thermal gas, cosmic rays and a magnetic field perpendicular to the gravitational 
accelerartion. Cosmic ray diffusion along the magnetic field is considered; cross 
field-line diffusion is supposed to be small and is ignored. We investigated two 
cases in our simulations. In the mechanical perturbation case we add a velocity 
perturbation parallel to the magnetic field lines, while in the explosional pertur- 
bation case we add cosmic ray energy into a sphere where the cosmic rays are 
injected. Linear analysis shows that the growth rate of the Parker instability 
becomes smaller if the coupling between the CR and the gas is stronger (i.e., 
the CR diffusion coefficient is smaller). Our MHD simulations of the mechanical 
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perturbation confirm this result. We show that the falling matter is impeded 
by the CR pressure gradient, this causes a decrease in the growth rate. In the 
explosional perturbation case, the growth of the magnetic loop is faster when the 
coupling is stronger in the early stage. However, in the later stage the behavior 
of the growth rate becomes similar to the mechanical perturbation case. 

Subject headings: cosmic rays — instabilities — ISM: magnetic fields — MHD 

1. INTRODUCTION 

It has been suggested that magnetic fields may play important roles for active phenom- 
ena in space, for example in astrophysical jets (e.g., Shibata & Uchida 1985; Matsumoto et 
al. 1996), solar activities (e.g., Priest 1982; Yokoyama & Shibata 2001), and active galaxies 
(e.g., Kuwabara et al. 2000). With such active phenomena, if a gas layer is supported by 
the horizontal magnetic fields against gravity, then the Parker instability may appear and 
can play an important role. Magnetic fields are also thought to play an important role in 
accretion disks (e.g., Stella & Rosner 1984; Kato & Horiuchi 1986) and galactic disks. For 
example, in galactic disks, the interstellar matter (ISM) is aggregated and grows to giant 
cloud complexes in spiral arms of galaxies via the Parker instability (e.g., Mouschovias 1974; 
Mouschovias et al. 1974; Elmegreen 1982a,b; Elmegreen & Elmegreen 1986). On the other 
hand, cosmic rays (CRs) may also play an essential role in the dynamics of the ISM, since 
it is recognized that the energy density of cosmic rays is of the same order as that of the 
magnetic field and turbulent gas motions (Parker 1969; Ginzburg & Ptuskin 1976; Ferriere 
2001). The importance of the effects of cosmic rays has been acknowledged, and a discus- 
sion concerning CRs was also given in the original work on the Parker instability (Parker 
1966). For studying the dynamics of CRs, there are several approaches. The particle-particle 
approach treats the plasma and the CRs as particles; the fluid-particle approach treats the 
plasma as a fluid and the CRs as particles; and the fluid-fluid approach treats the plasma and 
the CRs as fluids. The hydrodynamical approach to the Parker instability and Parker- Jeans 
instability without CRs has been done by nonlinear calculation (e.g., Matsumoto et al. 1988; 
Shibata et al. 1989; Chou et al. 2000; Kim et al. 2002). On the other hand, in spite of the 
suggestions of many astrophysical applications, there are very few papers on the evolution of 
the Parker instability with the effects of cosmic rays (Hanasz 1997; Hanasz & Lesch 2000). 
Hanasz & Lesch (2000) carried out calculations on the Parker instability induced by cos- 
mic ray injection from a supernova under the thin-flux-tube approximation, and suggested 
that the instability grows on shorter timescales, with the values of the diffusion coefficient 
decreasing. As the diffusion coefficient decreases, the diffusion speed of the CRs decreases. 
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Since the region where the CR energy is injected keep it for a long time, the dynamics is 
dominated by the interactions near the injection region. 

In this paper, we present the results of a linear analysis and MHD simulation of the 
Parker instability with the effects of cosmic rays, starting from an equilibrium two tempera- 
ture layered disk. We adopt the hydrodynamic approach for cosmic ray propagation (Drury 
& Volk 1981; Ko 1992). The paper is organized as follows. In §2, we present our physical 
assumptions and the equilibrium model. Linear stability analysis of the system is given in 
§ 3, and the numerical results are described in § 4. § 5 provides a summary and discussion. 



2. MODELS 



2.1. Assumptions and basic equations 

We investigate the Parker instability with the effect of the CRs in the galactic disk. The 
basic equations are the MHD equations combined with the CR energy equation, 
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where P g and P c are the gas pressure and the CR pressure, / is the unit tensor, 7 g and 7 C 
are the adiabatic index for the gas and the CRs, k\\ is the CR diffusion coefficient along the 
magnetic field, b is the unit vector of the magnetic field, f2 is the rotational angular frequency, 
and the other symbols have their usual meanings. In this model, self-gravity is ignored. The 
centrifugal force is assumed to be balanced by other forces (e.g., radial gravitational force 
of the galaxy). For simplicity, we ignore cross- field-line diffusion of the CRs. As a matter of 
fact, the ratio of perpendicular- to parallel-diffusion coefficient is quite small, say 0.02-0.04 
(Giacalone & Jokipii 1999; Ryu et al. 2003). We normalize these equations by the physical 
quantities related to the equilibrium model described in the next subsection. The unit of 
density and velocity are the density p an d sound speed C$o at the mid-plane of the galactic 
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disk in the equilibrium model. The unit of length is the scale height without magnetic field 
and CRs, H Q = Cg /(7 g <7 z ), and the unit of time is the sound crossing time over one scale 
height, Hq/Cso- The two dimensional calculation is carried out in the Cartesian coordinate 
system (x, z), where we adopted the approximation x = (p and z — z in the cylindrical 
coordinate system (r, tp, z) of the galactic disk, as did Mineshige et al. (1993) (see Figure 1). 
Moreover, the calculation is carried out only in the region over the mid-plane of the galactic 
disk. 



2.2. Equilibrium model 



We adopted the two temperature layered disk equilibrium model (Shibata et al. 1989) 
as the initial condition, 
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where the disk temperature is T = 10 4 K, the halo temperature is Th a i = 25 x 10 4 K, the 
height of the disk-halo interface is Zhaio = 900 pc, and the width of the transition layer is 
w tT = 30 pc. The magnetic fields are horizontal initially. The density, gas pressure, and CR 
pressure distributions are derived from the following equation: 
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subsequently, the total gas pressure scale height at z = (mid-plane of the galactic disk) 
is H = (1 + a + /3)Cg /(7g<7 2 ), where a, (3, and g z (> 0) are the initial ratio of magnetic 
pressure to gas pressure, the initial ratio of CR pressure to gas pressure, and the gravitational 
acceleration, respectively. In this paper, we only consider constant 7 g , 7 C , g z . In the following 
simulations, we pick 7 g = 1.05 and 7 C = 4/3, and set a — 1 and (3 = 1 initially. The system 
is initially homogeneous in the x-direction. For normalization, we take our units as follows: 
the unit of length is H = C| / (7 g g 2 ) = 50 pc, the unit of density is po = 1-6 x 10 -24 g cm~ 3 , 
the unit of velocity is C S o = 10 km s _1 , and the unit of time is H /C S o ~ 5 Myr, where the 
subscript denotes the value at the mid-plane of the galactic disk (z = 0). 



3. LINEAR STABILITY ANALYSIS 



3.1. Linearized equations 



We perform standard linear stability analysis on the set of equations (l)-(5). Since the 
unperturbed state depends on z only, the perturbed quantities are assumed to have the form 
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£ = (Sp, i SV X , i SV y , SV Z , SP g , SP C , SB X , SB y , —i SB Z ) and 

£ = £(z) exp(at + ik x + ik y ), 



(8) 



where a is the growth rate, k x and /c y are the wave numbers in x- and ^/-direction, respec- 
tively. For simplicity, instead of linearizing the energy equation, we assume an isothermal 
perturbation for the gas, 

5P g = C 2 s Sp, (9) 

where C s is the isothermal sound velocity. 

Although there are nine perturbed quantities, it turns out that there are seven algebraic 
relations among them, if we assume no cross-field-line diffusion of cosmic rays. First of all, 
the induction equation gives (SB X , SB y , 5B Z ) in terms of (5p,5V x ,5V y ,5V z ). The x- and y- 
momentum equations then give (SV X , 8V y ) in terms of (Sp, SV Z , SPt), where SP t = SP g + 5P c + 
B x SB x /4n. The nice consequence of no cross-field-line diffusion is that 5P C can be written in 
terms of (Sp,SV z ). (Note that the diffusion coefficient is also perturbed, because it depends 
on the direction of the magnetic field.) Thus p can be written in terms of (SV z ,SP t ). After 
some algebra, the continuity equation and the z-momentum equation become 
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We solve equation (10) to find the eigenmodes with given boundary values. Consequently, 
the problem converges to a boundary value problem, and the growth rate of the perturbation 
is obtained as an eigenvalue. 



3.2. Boundary conditions 

We assume that the disk is symmetric with respect to the mid-plane [z — 0). Under 
this assumption, the sign of 5V Z inverts beyond the mid-plane, on the other hand, the sign 
of SPt should not change. Hence, the perturbed value of y\ and 1/2 should be antisymmetric 
and symmetric about z — 0, respectively. On the other hand, the Matrix A in equation (10) 
nearly equals a constant in the region z ^> H and the WKB solution is applicable. Then, 
the asymptotic solutions, under the condition that y\ and 1/2 should vanish at large z, are 
written as follows (e.g., Horiuchi et al. 1988) 

y x = exp(Az), (19) 

(A - A n ) Vl = A 12 y 2 , (20) 

where 

A = - (A u + A 22 - y/(A n - A22) 2 + 4A 12 A 21 ) . (21) 

We solve the set of the two linear differential equations (eq.[10]) by the shooting method. We 
integrate equation (10) from the outer boundary by using the condition given by equation 
(19) and (20), to the inner boundary (z = 0) to obtain a trial value of a. This a is regarded 
as an eigenvalue when the value of y\ is small enough at z — 0. 



3.3. Result of linear stability analysis 



In this analysis, we take the value of the CR diffusion coefficient k\\, the ratio of CR 
pressure to gas pressure j3, and the rotational angular frequency Q as parameters. The value 
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of K|| is estimated as ~ 3 x 10 28 cm 2 s _1 (Berezinsky et al. 1990; Ptuskin 2001; Ryu et al. 
2003). In our units, 3 x 10 28 cm 2 s" 1 corresponds to 200. We thus take k\\ = 200, (3 — 1, and 
Q=0 as a fiducial case, and the ratio of magnetic pressure to gas pressure is taken as a = 1 
(initially). 

Figure 2a shows the dispersion relations for different k\\ values, a is the growth rate 
and k x is the wave number in the direction of the initial magnetic field. We set k y = for 
the linear analysis in this paper. In this calculation, we set (3 = 1, Q = 0. The growth rate 
becomes smaller as ku decreases from the fiducial value, ku = 200. This result well matches 
the one given by Ryu et al. (2003). The maximum growth rate <r max for different k\\ occurs at 
roughly the same k x . Figure 2b shows the dependence of the maximum growth rate cr max on 
«||- °max increases rapidly between < k\\ < 20, increases gradually between 20 < k\\ < 80, 
and finally, becomes almost constant when k\\ > 80. 

Figure 3a shows the dispersion relations for different (3, where (3 is the ratio of the CR 
pressure to the gas pressure. In this calculation, we set k\\ = 200 and Q = 0. The case of 
(3 = is identical with the Parker instability case without the effect of CRs. As (3 increases, 
the maximum growth rate <r max increases, and the wave number fc imax , at which the growth 
rate becomes maximum, increases as well. Moreover, <7 max and k xmax are roughly related 
linearly. Figure 3a also indicates that the short wavelength perturbations become unstable 
as (3 increases. Figure 36 shows the dependence of the maximum growth rate on f3. a max 
increases almost linearly in the range from f3 — to f3 — 1, and approaches a constant value 
as (3 increases beyond (3 — 1. 

Figure 4a shows the dispersion relations for different Q, where Q is the rotational angular 
frequency. In this calculation, we set Ku = 200. As Q increases, the growth rate of long 
wavelength perturbations decreases rapidly, and k x max becomes larger. Figure 46 shows the 
dependence of <7 max on (3 for different Q. 

4. TWO-DIMENSIONAL MHD SIMULATION 

4.1. Numerical Procedure and Boundary Conditions 

We solve the set of two-dimensional nonlinear, time- dependent, compressible ideal MHD 
equations, supplemented with the cosmic ray energy equation (eq.[l]-[5]), in Cartesian co- 
ordinates. We use the modified Lax-Wendroff scheme with artificial viscosity for the MHD 
part, and the Biconjugate gradients stabilized (BiCGStab) method for the diffusion part of 
the CR energy equation in the same manner as described in Yokoyama & Shibata (2001). 
The MHD code using the Lax-Wendroff scheme was originally developed by Shibata (1983) 
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and has been extended by Matsumoto et al. (1996) and Hayashi et al. (1996). To test the part 
of the code on diffusion, we did a test run on a simple CR diffusion problem (see appendix). 

We calculate only the region over the mid-plane of the galactic disk in the x-z plane, 
as shown in Figure 1. The x-direction corresponds to the azimuthal direction and the z- 
direction is in the direction of the rotational axis of the galactic disk. The partial derivative 
d/dy is neglected. We use the units in §2.2 for our simulations, and the equilibrium model 
in §2.2 is used as the initial equilibrium background. We then apply two types of pertur- 
bation, mechanical perturbations and explosional perturbations. In the case of mechanical 
perturbations, we add small velocity perturbations of the form 

(22) 

within the finite rectangular region, < x < A/2 and AH < z < 8H , where A = 20 H 
is the horizontal wavelength of the small velocity perturbation, which is nearly equal to 
the most unstable wave length derived from the linear analysis in the km = 200 case. The 
size of the simulation box is (x max x z max ) = (80H x 187 H ), the number of grid points is 
(N x , N z ) = (101, 401), the grid size is Ax = 0.8H and Az = 0.15H , when < z < 25.0# , 
and increases with z otherwise. We assume a symmetric boundary for x = and z — 0, 
and a free boundary for x = x max and z = z mSuX . On the other hand, in the explosional 
perturbation case, the cosmic ray energy (~ 10 50 ergs) is put into a cylindrical region which 
is placed at (x, z) = (0, 6H ) with volume V cxp = vrr^ xp y exp , where r exp = 25 pc = 0.5H , 
y CX p = 50 pc = H . The size of the simulation box is (rc max x z max ) = (90H x 187iJ ), 
the number of grid points is (N x , N z ) = (301, 401), the grid size is Ax = Az = 0.15if , 
when < x < 35Hq and < z < 25Hq, and increases with x and z otherwise. We assume 
symmetric boundaries for x = and z — 0, and a free boundary for x = x max and z = z max . 



v x = 0.05 x Cgo sin 
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4.2. Numerical Results for The Mechanical Perturbation Case 

In this subsection, we show the results for the mechanical perturbation case. In order to 
examine the effect of the cosmic ray diffusion coefficient on the instability, we consider three 
different diffusion coefficients: k\\ = 10, k\\ = 40, and k\\ = 200. As can be seen from Figure 
2, km = 10, K\\ = 40, and km = 200 correspond to a small, medium, and high growth rate, 
although the difference of the growth rate corresponding to k\\ = 40 and that to k\\ = 200 is 
small. 

Figure 5a-c show the time evolution of the CR pressure distribution, the magnetic field 
lines, and the velocity vectors for the of Ky = 200, Ky = 40, and Ky = 10 models. The gray 
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scale contour shows the CR pressure distribution, white curves show the magnetic field lines, 
and white vectors show the velocity vectors, where a white arrow at the upper right corner 
shows a reference velocity vector (= 5 x Cso)- The middle column shows the result of the 
three models (ku = 200, ku = 40, and ku = 10) at the end of their linear growth (i.e., t = 28, 
t = 30, t = 36, respectively). These times are decided from Figure 6, which shows the time 
evolution of V x in each model. The initial perturbations grow to form the characteristic loop 
like structures. In the linear phase, the form of the magnetic loop is almost the same in each 
model. In contrast, in the non-linear phase, its form can be very different. In the k\\ = 200 
model, the shape of the magnetic loop is a beautiful omega shape. In the k\\ = 40 model, 
the omega shape is distorted in the outer region and a double-loop forms in the inner region. 
In the K|| = 10 model, the major feature is the double-loop. 

In order to compare the results of the MHD simulation and the linear analysis, we 
examine the temporal variation of V x at a particular point. Figure 6 shows the time evolution 
of V x at (x, z)—(5.6, 7.95). The solid curves show the growth rate of V x (normalized to the 
sound velocity C s q) in each model: ku = 200, k» = 40, and k» = 10. The broken lines, 
LI for k\\ = 200, L2 for k\\ = 40, and L3 for k\\ = 10, show the growth rate obtained from 
the linear analysis. The dash-dotted line shows the initial Alfven speed. The inclination 
obtained from the linear analysis agrees well with that obtained from the MHD simulation. 
The speed reaches the Alfven speed and is saturated in the k\\ = 200 model. On the other 
hand, the speed is saturated below the Alfven speed in the k\\ = 40 and k\\ = 10 models, 
and the saturated speed decreases as k\\ decreases. 

In Figure 7, the three panels in the first row show the CR pressure distribution (gray 
scale contour), velocity vectors (white arrows), and a magnetic field line (white curve), at 
the end of the linear growth phase for the three models: (a) k\\ = 200 at t — 28t , (b) k\\ = 40 
at t — 30t , (c) K|| = 10 at t — 3Gt . The arrow at the upper-right corner shows a reference 
velocity vector equal to 5 times a unit velocity vector. As the value of k\\ decreases, the 
expansion speed of magnetic loop becomes slower. The three panels in the second row show 
the CR pressure values along a magnetic field line depicted in the first row. In the ku = 200 
model, the CR pressure distribution in the magnetic loop (x < 15) is nearly uniform. As 
K|| decreases, the CR pressure profile develops some structures. In addition, as k\\ decreases, 
the CR pressure decreases at the top of the loop (x m 0) but increases at the foot point of 
the loop. The horizontal axis l U is the distance along a magnetic field depicted in the first 
row. L = 64.5 (k\\ = 200), L = 61.1 (k\\ = 40), and L = 57.8 (k\\ = 10) at x = 50. The 
three panels in the third row show the density distributions, log 10 (p/p ), along a magnetic 
field line depicted in the first row. The distributions are very similar in each model, but we 
can recognize the difference of the distribution in the magnetic loop part. In the k\\ = 200 
model, the density in the loop is lower than that in the others. In the other models, the 
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density increases in the loop as k\\ decreases. The three panels in the fourth row show the 
absolute value of the velocity at the position of the magnetic field line depicted in the first 
row. The absolute value of velocity is large where the curvature of the magnetic loop is large 
(except at the foot point). Its maximum value becomes smaller as k» decreases. The three 
panels in the fifth row show the velocity component along the magnetic field line depicted 
in the first row. When the value is positive /negative the velocity is parallel/anti-parallel to 
the magnetic field. The falling speed of matter along a magnetic field line is supersonic near 
the foot point of the magnetic loop in the rey = 200 and k\\ = 40 models. On the other hand, 
it is subsonic in the k\\ = 10 model. 

Figure 8a-c show the time evolution of V z (velocity along the z-axis) sliced at x = for 
the k\\ = 200, K|| = 40, and k\\ = 10 models. In each model, V z becomes maximum at the 
position of the magnetic loop near the end of the linear growth, t ~ 29 for k\\ = 200, t ~ 31 
for k\\ = 40, and t ~ 36 for k\\ = 10. Subsequently, V z becomes large in the halo region, 
because the halo is pushed upward by the growing magnetic loops. The growth speed of 
the magnetic loop finally falls, and the growth is impeded similarly to the case without the 
effect of CRs (Kato, Fukue, & Mineshige 1998, chap. 17). 

4.3. Numerical Results for The Explosional Perturbation Case 

In this subsection, we show the results for the explosional perturbation case. In this 
case, we set a high CR pressure region at (x, z)—(0, 6) with radius 0.5H as our initial 
perturbation. 

Figure 9 shows the time evolution of the CR pressure distribution (gray scale contour), 
the magnetic fields (white curves), and the velocity vectors (white arrows) for the models 
(a) km = 10, and (b) k» = 80. It is recognized that the growth speed of the magnetic loop in 
the k\\ = 10 model is faster than the k» = 80 model at t — 12, contrary to the result for the 
mechanical perturbation case. Subsequently, at t — 24, the growth of the magnetic loop in 
the AC|| = 80 model overtakes the k» = 10 model. The expansion speed of the magnetic loop 
becomes very slow in the k\\ = 10 model, while it is still very fast in the k\\ = 80 model. 

In Figure 10, the two panels in the first row show the initial CR pressure distribution 
(gray scale contour), an initial magnetic field line (white line) in the models k» = 10 (the 
left column) and k\\ = 80 (the right column). The two panels in the second row show 
the time evolution of the CR pressure distribution along a magnetic field line depicted in 
the first row. The field line threads through the explosion region. The dotted line shows 
the initial distribution of the CR pressure, which is highly localized near x — 0. In the 
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K|| = 10 model, the high CR pressure region is localized for a relatively long time, because 
the diffusion speed is slow. The initially localized high CR pressure pushes the matter in the 
x-direction along the magnetic field lines rather effectively in the case of strong coupling (i.e., 
small diffusion coefficient). Thus the density drops rather rapidly, while the CR pressure 
decreases slowly in the initial phase (t < 2.5). When the magnetic loop penetrates into the 
low CR pressure region (t > 2.5), the CR pressure inside the loop diminishes faster, partly 
because the magnetic tube has expanded and partly because the cosmic ray is carried by the 
downward flow of the matter. As matter accumulates at the foot point of the magnetic loop, 
the magnetic tube becomes thinner, and the CR pressure builds up, because of the small 
When a significant CR pressure gradient has been built up against the infall of matter, 
the growth of the instability slows down. This can be confirmed by the small differences of 
the density and the CR pressure distributions between t — 20 and t = 24. In the k» = 80 
model, diffusion is more important. The high CR pressure region disappears as cosmic rays 
diffuse along the magnetic field. At t — 15, the CR pressure is rather uniform throughout 
the whole region. The CR pressure gradient between the top of the loop and the foot point is 
significant only after t = 24, and the growth of the instability will not slow down until then. 
The two panels in the third row show the time evolution of the density distribution along 
a magnetic field line depicted in the first row. The field line threads through the explosion 
region. In the k\\ = 10 model, matter is drained rapidly by the CR pressure gradient and a 
large drop in density occurs near the top of the magnetic loop in a short time (at t = 2.5). 
As time proceeds, the matter accumulates at the foot point of the loop where a high density 
region is formed. The draining rate of matter near the top of the loop reduces. In the 
K|| = 80 model, the density near the top of the magnetic loop decreases very slowly, until 
about t = 15. After that, the draining rate accelerates, and the density near the top of the 
loop becomes smaller than that of the rey = 10 model at t — 24. The two panels in the last 
row show the CR pressure distribution, the velocity vectors, and the magnetic field line at 
t = 24. We should point out that, in order to emphasize the CR pressure distribution near 
the depicted magnetic field line, the gray scale used in Figure 10 is different from that used 
in Figure 9. 



5. SUMMARY AND DISCUSSION 

Using linear analysis and a time-dependent non-linear calculation, we studied the Parker 
instability (or magnetic buoyancy instability) with the effect of cosmic rays. Several works on 
the linear analysis of the Parker instability with the effect of CRs have been published (e.g. 
Hanasz & Lesch 1997; Ryu et al. 2003). In Hanasz (1997) the CR energy equation (including 
diffusion) was not solved. In Ryu et al. (2003) the effect of rotation is not included, and only 
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two cases of CR pressure were described (one was without CR (3 — 0, and the other was with 
the same unperturbed CR and gas pressures (3=1). Since Ryu et al. (2003) showed that 
the effect of cross-field-line diffusion of CRs is negligible in the context of ISM, we neglected 
the effect of cross-field-line diffusion in our analysis for simplicity. 

In the linear analysis, the growth rate becomes larger as the CR diffusion coefficient 
along the field line k\\ increases, and is saturated at large k\\. This result is consistent with 
the result by Ryu et al. (2003). The growth rate also becomes larger when the initial ratio of 
the CR pressure to gas pressure (3 increases, and is saturated at large (3. This is consistent 
with Ryu et al. (2003), except for some slight differences. In our results, the maximum 
growth rate of the normal Parker case ((3 = 0) is almost half of the fiducial case ((3 — 1), 
and the critical wave number, over which the instability is stabilized, of the normal Parker 
case is about 0.7 times of the fiducial case. In Ryu et al. (2003), the maximum growth rate 
of the normal Parker case is less than half of the (3 = 1 case, and the critical wave number 
of the normal Parker case is about half of the (3 = 1 case. The differences perhaps come 
from how the normalization was taken. In fact, we succeeded in producing their results by 
taking the same scale height under the same equilibrium condition. The scale height of the 
disk, H = (1 + a + /3)C s 2 o/(7gfi'z), changes with the values of a and (3, and it takes the value 
H = 2Cg / '(7 g <7z) in the normal Parker case (a — 1, (3 — 0), and H = 3Cg /(7 g g 2 ) in the 
fiducial case (a — (3 — 1). We allowed the change of the scale height because we preferred 
not to change the gravitational acceleration. This is the reason why we took the unit of 
length as H = C^/(^ g g z ). The effect of the rotation stabilizing the Parker instability is 
similar to the case without CRs. The k xmax increases as the increases. Our result for the 
effect of rotation is consistent with that by Hanasz & Lesch (1997), except for the difference 
in the region of small wave number. The growth rate with rotation is small near k x = 0. It 
increases linearly with the wave number in Hanasz & Lesch (1997). However, in our result 
the growth rate increases faster than linear, and this is also observed in the normal Parker 
case (see Kato, Fukue, & Mineshige 1998, chap. 17). 

In the MHD simulation, we showed that the growth rate of the instability becomes 
smaller when the diffusion coefficient k» becomes smaller, which agrees well with the result 
of the linear analysis. At the end of the linear growth, the morphology of the magnetic loop 
developed from the initial perturbation is more or less the same in the three models (/cy = 200, 
k\\ = 40, and k\\ = 10) studied in the mechanical perturbation case. However, in the non- 
linear stage, the magnetic loop in different models develops into different morphologies. 
From the distribution of CR pressure, density, and velocity along a magnetic field line at 
the end of linear growth, we found several characteristics. In the case of small diffusion 
coefficient (e.g., k\\ = 10, i.e., the coupling between the CRs and the gas is strong), the CR 
pressure distribution is rather non-uniform. CRs tend to accumulate near the foot point of 
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the magnetic loop, and the CR pressure gradient force towards the top of the loop becomes 
larger. The falling motion of matter is then impeded by the CR pressure gradient force, 
and the growth rate of the Parker instability decreases. In the case of a large diffusion 
coefficient (e.g., k» = 200, ku = 40), the falling speed of matter along a magnetic field 
line exceeds the speed of sound, and a shock is formed near the foot point of the magnetic 
loop. Moreover, the CR pressure distribution along a magnetic field line in the cases of large 
diffusion coefficient (e.g., rey = 200) reminds us of the profile of CR pressure in cosmic-ray- 
modified shocks (e.g., Drury & Volk 1981; Ko et al. 1997). The linear growth rate in the 
simulations agrees well with that in the linear analysis. We also found that the speed along 
the disk is saturated at the initial Alfven speed. This result agrees with that in the normal 
Parker instability (i.e., without CRs, Matsumoto et al. 1988). The unperturbed state has the 
scale height H = (C s 2 +/3Cg + V£/2)/g z . When the Parker instability takes place, the gas falls 
down along the magnetic field lines. The CR pressure tends to distribute uniformly along a 
magnetic field line (at least in the case of large diffusion coefficient), and its contribution to 
the scale height disappears. Thus the scale height along the magnetic field lines settles to 
H' = Cg/g z at later times. The released gravitational energy in the form of kinetic energy 
per unit mass is estimated as V^/2. Hence, we obtain the same results as the normal Parker 
case, even when the effect of CRs is included. 

The explosional perturbation case has been studied by Hanasz & Lesch (2000). They 
stated that the smaller the diffusion coefficient the larger the growth rate of the instability. 
This trend is the opposite of what we found from linear analysis and simulation in the 
mechanical perturbation case. We thus computed the explosional perturbation case for 
a longer time. Our result showed that the growth rate is larger in the smaller diffusion 
coefficient model only in the early stage. The growth rate becomes smaller when compare 
with that of the large diffusion coefficient model in the later stage. The growth of instability 
is suspended by the CR pressure gradient force interfering with the falling motion of the 
matter in the small k\\ model, while the magnetic loop can grow up to larger scales in the 
large k\\ model. 

Numerical computations were carried out on VPP5000 at NAOJ. TK and CMK are 
supported in part by the National Science Council, Taiwan, R.O.C., under the grants NSC- 
91-2112-M-008-006, NSC-90-2112-M-008-020, NSC-91-2112-M-008-050. 
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A. TEST CALCULATION 

In this appendix we show the result of a simple CR diffusion problem to test the diffusion 
part of our numerical code. We solved the following diffusion equation (i.e., cosmic ray energy 
equation [4] with V = 0), 

-^ = V.[«||66.VS C ], (AI) 

where E c is the CR energy, and ku is the diffusion coefficient along the magnetic field. We 
considered a uniform magnetic field with x-component only, i.e., B = (B x ,0). The test 
calculation itself is two dimensional in the x-z plane but the content is the same as a one 
dimensional calculation in the x-direction, because we just considered constant ku and the 
initial condition depended on x only. We took the same initial condition as that used in 
Hanasz & Lesch (2003), 



Eqo = A exp 



,y»2 



(A2) 



where w = V50 is the initial half width of the Gausian profile, x p is the distance from 
the central point of calculation region, and A = 10 is the value at x p = 0. The number of 
grid points used in this calculation is (JV X , iV 2 ) = (400, 100), and k\\ = 100. Figure 11 shows 
the initial distribution and the distribution at t — 9.4. In the right panel (t = 9.4), the 
solid curve shows the analytical solution and the squares show the numerical solution. The 
numerical solution completely matches the analytical solution. 
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Fig. 1. — Schematic picture of the simulation model and simulation box. 

Fig. 2. — (a) Dispersion relation for the Parker instability with the effect of CRs at different 
K\\. a is the growth rate of perturbation and k x is the wave number along the direction of 
the magnetic field in the unperturbed state, (b) The dependence of the maximum growth 
rate <r max on k\\. 

Fig. 3. — (a) Dispersion relation for the Parker instability with the effect of CRs at different 
(3, where (3 is the initial ratio of CR pressure to gas pressure, (b) The dependence of the 
maximum growth rate cr max on (3. 

Fig. 4. — (a) Dispersion relation for the Parker instability with the effect of CRs at different 
Q, where fl is the rotational angular velocity of the disk. (b)The dependence of the maximum 
growth rate cr max on (3 for different Q. 

Fig. 5. — Time evolution of the CR pressure distribution, magnetic field lines, and velocity 
vectors in the models (a) km = 200, (b) km = 40, and (c) km = 10. The gray scale, the white 
curves, and the white vectors show the CR pressure distribution, magnetic field lines, and 
velocity vectors. The units of length and time are 50pc and 10 6 years, respectively. 

Fig. 6. — Comparison of growth rate between the results of linear analysis and the results 
of MHD simulation. The broken lines, LI, L2, and L3 show the power law relation for the 
models k\\ = 200, Ky = 40, and Ky = 10. The dash-dotted line shows the initial Alfven speed. 

Fig. 7. — The CR pressure, the density, the absolute of velocity and the speed along a 
magnetic field line depicted in the first row at the end of linear growth in the models (a) 
Ky = 200, (b) k\\ = 40, and (c) k\\ = 10. L is the distance along the magnetic field line. 

Fig. 8. — Time evolution of V z sliced at x = in the models (a) Ky = 200, (b) Ky = 40, and 
(c) Ky = 10. 

Fig. 9. — Time evolution of the CR pressure distribution (gray scale), magnetic fields (white 
curves), and velocity vectors (white arrows) in the case of explosional perturbation. 

Fig. 10. — Time evolution of the CR pressure distribution, and the density distribution along 
a magnetic field initially threading the region where the CR energy was injected. The left 
column is km = 10, and the right column is km = 80. 
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Fig. 11. — Test result of a simple CR diffusion problem. 
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